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ABSTRACT 

Path  and  signal  processing  corrections  made  to  amplitudes  give  magnitudes  mb  and  Ms,  that  in  principle,  have  only 
source  information.  We  propose  the  addition  of  a  model-based  magnitude  prediction  as  a  correction  term,  that  is,  a 
source  correction  in  addition  to  path  and  signal  processing  corrections  under  the  null  (HO)  hypothesis  that  a  seismic 
event  is  a  single-point  fully  contained  explosion.  This  additional  correction  removes  explosion  source  information 
from  the  Ms  versus  mb  discriminant  with  the  remaining  source  information  represented  as  a  constant.  There  are 
effects  such  as  depth,  focal  mechanism,  and  local  material  properties  that  cannot  easily  be  determined  and 
mathematically  included  in  amplitude  corrections.  We  develop  a  mathematical  model  to  capture  these  near  source 
effects  as  random  (unknown)  giving  an  error  partition  of  two  sources:  model  inadequacy  and  station  noise.  This 
mathematical  model  is  the  basis  for  two  new  Ms  versus  mb  discriminant  formulations.  Both  methods  are  designed  to 
utilize  advances  in  source  physics  theory  by  using  new  source  models  to  predict  (correct)  observed  surface  wave 
magnitudes.  The  network  average  discriminant  formulation  includes  source  model  error  and  magnitude  correlation 
in  the  standard  error  (SE)  of  the  discriminant,  effectively  placing  a  lower  bound  on  the  discriminant  SE  and 
accounting  for  scaling  between  Ms  and  mb.  This  new  property1  of  the  Ms  versus  mb  discriminant  correctly  reduces 
the  SE  only  through  network  averaging  of  magnitudes.  A  second  discriminant  formulation  derived from  order 
statistics  theory  is  potentially  robust  to  station  coverage  for  an  event  by  basing  discrimination  on  the  maximum 
observed  station  Ms  value 
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OBJECTIVES 

When  completed,  this  research  will  leverage  advances  in  Ms  signal  processing  into  significant  reformulations  of  the 
Ms  versus  mb  discriminant  to  include:  1 )  a  new  correction  for  source  in  addition  to  established  path  corrections 
under  the  null  hypothesis  that  a  seismic  event  is  an  explosion,  2)  a  new  sources-of-error  model  that  correctly 
partitions  total  error  between  two  components  -  model  and  background  noise,  and  provides  the  correct  formulation 
for  scaling  between  Ms  and  mb,  3)  and  a  new  Ms  versus  mb  discriminant  formulation  that  can  be  robust  to  surface 
wave  radiation  pattern  (station  coverage)  for  earthquakes.  We  review  progress  to  date  on  this  research. 

RESEARCH  ACCOMPLISHED 

An  instrument-corrected  amplitude  spectrum  can  be  thought  of  as  a  convolution  between  the  source  and  the  path.  In 
the  frequency  domain  this  can  be  mathematically  represented  as 

A(c,A)  ,WA)P(.)0(,,A)  j  ()) 

where  S  is  the  source  spectrum,  G  is  geometrical  spreading,  P  is  the  frequency-dependent  site  effect,  and  B  is  the  an- 
elastic  attenuation  with  function  arguments  epicentral  distance  A  and  angular  frequency  w.  I  lere  we  have  split  the 
path  effect  into  three  components:  1 )  a  frequency  independent  geometrical  spreading  component,  2)  a  range 
independent  and  frequency  dependent  site  effect,  and  3)  an  an-elastic  attenuation  component. 

The  logarithm  of  both  sides  of  Equation  (1)  gives 

logA(^,A)  logS(cj)  +  log  (7(A)  |  log/V)  1  log  13(a),  A)  (2) 

In  general,  to  remove  source,  size  and  path  trends  in  the  data,  we  can  correct  the  observed  log  ^(/co.A)  spectrum  so 
that 


log  Ac(u;.  A)  log,4„(uJ.  A)  -  log  ,4(u;,  A^  ^ 

where  log  /!,{co,A)  is  the  corrected  spectrum.  Equation  (3)  gives  corrected  magnitudes  (residuals)  Y  that  are  then 
used  to  construct  discriminants.  Specifically,  from  Equation  (3),  the  corrected  magnitude  Y  is  a  log  observed 
magnitude  minus  a  model  prediction  of  magnitude  under  the  hypothesis  that  the  event  was  generated  by  source  S(to) 
(e.g,  single-point  fully  contained  explosion).  The  residuals  Y provide  data  (evidence)  to  test  this  hypothesis. 
Parseval's  theorem  provides  the  signal  processing  for  magnitude  calculations  (total  phase  energy)  in  the  time 
domain. 


A  baseline  source  model  5(co)  will  forecast  estimates  of  mb  and  Ms  for  an  event  in  question.  Fully  developed,  the 
source  model  S( co)  is  a  linear  superposition  of  two  axisymmetric  force  systems  for  explosions  conducted  in  hard 
rock.  First,  a  monopole  represents  the  explosion  as  a  spherical  source.  Second,  a  compensated  linear  vector  dipole 
(CLVD)  with  vertical  axis  in  extension  represents  shock-induced,  deep-seated  tensile  failure,  a  form  of  source 
medium  damage  (see  Patton  and  Taylor  (2008),  and  Ben-Zion  and  Ampuero  (2009)).  Both  sources  are  assumed  to 
be  time  and  space  coincident,  and  have  a  common  Mueller  and  Murphy  (1971)  source  function.  The  monopole 
source  alone  is  used  to  forecast  mb.  Ms  is  forecast  by  half-space  calculations  of  Rayleigh  wave  interference  between 
the  monopole  and  CLVD  force  systems  (Patton  and  Taylor  (2008)). 

Estimates  of  mb  and  Ms  and  error  bounds  are  derived  from  statistical  analysis  of  model  parameters.  With  density 
and  P  velocity  fixed  for  hard  rock,  the  model  has  four  parameters:  Poisson  ratio  v,  depth  of  burial  h,  yield  IF,  and  a 
low-frequency  CLVD  source  strength  to  monopole  parameter  K  developed  in  Patton  and  Taylor  (2008).  Independent 
observations  constrain  or  bracket  values  for  each  parameter.  Such  observations  might  be  drawn  from  site  geology 
and  tectonic  history  (v),  previous  testing  history  (mature  versus  new  test  sites;  h,  K),  and  seismic  measurements, 
both  regional  and  teleseismic,  of  source  size  (IF).  The  two  discriminants  developed  here  utilize  a  prediction  of  mb 
and  Ms  derived  from  a  reasonable  suite  of  models  covering  the  parameter  space.  Thus,  the  discrimination  methods 
developed  in  this  paper  have  the  potential  to  provide  monitoring  capability  in  settings  with  no  explosion  calibration 
data. 
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Consistent  with  the  established  approach  to  Ms  versus  mb  discrimination,  we  develop  the  mathematics  of  two  new 
discriminant  formulations  -  a  multi-station  discriminant  constructed  from  the  network  average  of  the  corrected 
magnitudes  F,  and  a  multi-station  discriminant  constructed  from  the  maximum  of  station  Ms  values.  Both 
discriminants  are  built  from  random  effects  analysis  of  variance  (Scheffe  (1959),  Searle  (1971)  and  Searle  et  al. 
(1992))  which  has  been  applied  to  other  path  correction  theories  in  seismology  (e.g.,  Chen  and  Tsai  (2002),  Tsai  and 
Chen  (2003)  and  Tsai  et  al.  (2006)).  We  model  any  remaining  physical  structure  in  corrected  magnitudes  as  a  source 
bias  plus  two  random  effect  components  -  model  inadequacy  and  station  noise.  This  approach  to  Ms  versus  mb 
discriminant  formulation  properly  forms  the  discriminant  standard  error  with  these  two  variance  components.  Model 
inadequacy  decreases  with  scientific  advances  in  source  and  path  correction  theory  and  improved  calibration,  and 
consistent  with  signal  processing  theory,  station  noise  is  reduced  through  station  averaging  in  the  network  average 
formulation.  The  conceptual  motivation  for  the  new  discrimination  formulations  is  illustrated  in  Figure  I. 
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(c)  Ms  Magnitudes  Corrected  f»r  Y icld Magnitude  and  Smuree 

Figure  1.  Conceptual  illustration  of  corrections  applied  to  Ms  magnitudes,  conditional  on  known 
Yield  Magnitude.  Normal  curves  represent  the  distribution  of  station  Ms  values. 

Figure  1(a)  illustrates  uncorrccted  Ms  magnitudes.  Figure  1(b)  illustrates  the  event 
populations  after  correction  for  Yield  Magnitude  giving  a  significant  reduction  in  model 
error.  Figure  1(e)  illustrates  the  event  populations  after  correction  for  Yield  Magnitude 
and  additionally  a  Source  correction  (observed  minus  prediction).  Note  the  additional 
reduction  in  model  error  and  the  significant  source  population  separation. 


364 


2010  Monitoring  Research  Review:  Ground-Based  Nuclear  Explosion  Monitoring  Technologies 


Discriminant  Formulation  with  Model  Inadequacy  and  Station  Noise  Error  Partition 

Established  signal  processing  treats  amplitudes  as  lognormal  distributed  and  therefore  magnitudes  are  normally 
distributed.  The  conceptual  representation  of  the  proposed  model  is 

Y  log(oorrc‘cted  amplitude)  f’(soura')  I  Event  I  Noise  ^4) 

where  C(source)  is  a  source  constant,  Event  is  a  zero  mean  random  effect  that  varies  from  event  to  event  and 
represents  model  inadequacy  from  effects  such  as  depth,  focal  mechanism,  local  material  property  and  apparent 
stress  variability,  and  Noise  represents  measurement  and  ambient  noise,  also  with  zero  mean.  This  approach  results 
in  C(source)  terms  for  earthquakes  and  explosions,  with  large  differences  in  these  terms  indicative  of  good 
discrimination  ability.  In  subsequent  development,  the  hypotheses  test  that  an  event  is  a  single  point  fully  contained 
explosion  will  be  mathematically  represented  in  terms  of  these  C(source)  parameters.  Equation  (4)  implies  the 
expected  value  of  Y  is  E  { F)  =  C(source).  In  the  development  it  is  assumed  that  good  signal  processing  practices 
have  been  applied  giving  high  confidence  in  the  quality  of  the  observed  amplitudes  and  calculated  magnitudes. 

For  the  mathematical  statistics  formulation  of  Equation  (4),  define  the  random  variable  to  be  a  magnitude 
residual  for  source  /=  0, 1  (explosion,  earthquake),  event  j  and  station  k  (observed  data  are  denoted  v,y*).  The  statistical 
linear  model  representation  of  Equation  (4)  is  then 

^  ijk  Hi  4  Ej  1  Hij)k  j  1,2,  ...rrij  k  1,2,  ...ri;;-  ^ 


Analogous  to  Equation  (4),  Equation  (5)  reads  V  equals  a  source  constant  p,  plus  a  random  event  adjustment  £) 
(source  and  path  model  inadequacy)  plus  a  station  noise  adjustment  e^*.  The  subscript  notation  (ij)k  for  ef,yp 
specifies  that  observed  station  noise  (k)  is  different  for  each  source  and  event  combination  (ij).  Equation  (5)  is  a 
standard  mixed  effects  (random  and  fixed)  linear  model. 

The  Et  are  modeled  as  independent  normal  random  variables  with  zero  mean  and  variance  f .  The  c^x  are 
independent  normal  random  variables  with  zero  mean  and  variance  cr.  £)  and  c(,j j*  are  independent  across  all 
subscripts.  This  assumption  is  consistent  with  effects  local  to  the  source  being  uncorrelated  with  station  noise. 

Statistical  Properties  of  a  Network  Magnitude  Residual 

For  source  /  and  event  j,  denote  the  1  x  n,j  vector  of  magnitude  residuals  for  n,j  stations  as  Y',r  Then  Y',j  is 
multivariate  normal  with  1  x  iij,  mean  vector  0',y=(p/,  p,,...,  p,j  and  n,j  x  /i/;  covariance  matrix 


/r2  +  rr2 

r2 

r2 

T~  \ 

r2 

r2  I  a2 

T2 

T2 
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t2  (  a2  t2 

V  r2 
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T-  T2  ^  o'2 J 

The  network  magnitude  residual  is  Y,=  HYjjln,j  is  normal  with  mean  p,  and  standard  error  r+crV/j^.  Omitting  the 
term  Ej  in  Equation  (5)  implies  that  the  magnitude  residual  at  a  station  is  p,  plus  station  noise.  This  model 
formulation  is  fundamentally  inconsistent  with  seismic  observation.  The  standard  error  of  Y,j  with  Ej  removed  from 
Equation  (5)  is  cf/n, ,■  ( r=0)  and  decreases  as  the  number  of  stations  observing  an  event  increases.  By  not 
including  the  term  £,,  effects  such  as  depth,  focal  mechanism,  local  material  property  and  apparent  stress  variability 
are  not  accounted  for  in  the  theoretical  model  representation  of  a  magnitude,  and  this  bias  can  never  be  diminished 
with  a  network  average  calculation.  Equation  (5)  captures  source  and  path  model  deficiencies  as  a  random  effect  (Ej) 
and  compensates  for  them  as  a  component  in  the  standard  error  of  a  discriminant.  Also,  the  lower  bound  of  a 
magnitude  standard  error,  derived  from  Equation  (5),  is  non-zero  and  therefore  consistent  with  realistic  seismic 
observation. 

Another  important  property  of  this  model  is  that  a  magnitude  residual  for  a  single  event  is  correlated  across  stations. 
This  correlation,  r/(r+cr)  implies  that  large  adjustment  Ej  increases  correlation  between  stations  because  this 
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random  adjustment  is  applied  to  all  stations  observing  an  event,  that  is,  the  station  residuals  stochastically  move 
together.  Small  adjustment  £)  implies  the  correction  model  is  good  and  is  conceptually  equivalent  to  stations  with 
incoherent  noise.  Small  adjustment  £j  also  implies  r  is  small  and  the  standard  error  of  T,;  is  fundamentally  reduced 
through  network  averaging. 

The  corrected  magnitudes  mb  and  Ms  can  be  correlated  if  they  are  both  corrected  for  size  in  the  source  spectrum  5 
with,  for  example,  Mr.  If  Ms  is  corrected  for  size  with  mb,  then  these  magnitudes  will  be  uncorrelated.  From 
Equation  (5)  the  random  mechanism  causing  this  correlation  is  modeled  with  the  £}  terms  between  magnitudes 
because  these  errors  perturb  the  signal  equally  at  all  stations  observing  an  event,  and  also  the  measured  amplitudes. 
This  correlation  is  the  mathematical  mechanism  that  captures  any  scaling  between  magnitudes.  Extending  Equation 
(5)  to  two  corrected  magnitudes  Y I  yt  (=  Ms,^)  and  Y2yk  (=  mbyt)  with  the  Eij  and  E2t  terms  between  magnitudes 
correlated  (p)  provides  the  statistical  model  to  calculate  the  standard  error  of  the  multi-station  discriminant  given  in 
the  section  that  follows. 

Network  Average  Discriminant  Formulation 

The  Ms  versus  mb  discriminant  is  constructed  from  a  bivariate  normal  model  of  network  magnitude  residuals  Ms 
and  mb,  with  bias  constants  Pag/  and  pm*/,  ;=0,I  (explosion,  earthquake).  For  an  explosion,  the  network  magnitude 
residual  Ms  is  normal  with  mean  piWjfl  and  standard  error  f  Kid  &  mJ> Ufa  and  mb  is  normal  with  mean  pm*0  and 
standard  error  f  mtdcf  mb!nmb-  With  the  inclusion  of  correlation  between  model  errors  Esuj  and  Emy,  the  standard 
error  of  the  Ms  versus  mb  discriminant  is  then 
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(7) 


for  both  earthquakes  and  explosions.  Equality  of  the  standard  error  for  both  source  types  is  an  important  model 
property  because  unlike  discrimination  analysis  with  unequal  source  variability,  it  ensures  that  an  earthquake  with 
an  unusually  strong  earthquake-identifying  discriminant  value  will  not  be  identified  as  an  explosion.  This  can  occur 
with  discrimination  methods  that  model  source  variance  as  unequal  -  explosion  calibration  data  can  exhibit 
variability  that  is  significantly  larger  than  earthquake  calibration  data. 

In  an  operational  setting,  physical  correction  theory  will  never  be  able  to  adjust  amplitudes  for  all  local  systematic 
effects.  As  discussed  previously,  these  local  systematic  effects  are  modeled  as  random,  moving  out  to  all  stations 
(and  therefore  amplitudes),  hence  the  correlation  between  a  network  of  stations.  As  physical  (source  and  path) 
corrections  improve,  the  model  inadequacy  terms  Ej  will  be  small,  giving  small  values  of  r,  and  the  covariance 
between  amplitudes  will  be  small.  In  the  limit,  this  conceptually  gives  station  discriminant  scatter  plots  for 
explosions  and  earthquakes  (populations)  that  are  small  shotgun  patterns  of  data. 

Centering  relative  to  the  explosion  population  means,  the  standardized  discriminant  is 
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which  is  centered  at  zero  for  explosions  and  has  a  non-zero  center  for  earthquakes.  The  advantage  to  centering 
relative  to  explosions  is  consistency  with  the  monitoring  position  to  assume  all  events  are  explosions  and  then  prove 
otherwise  with  seismic  signatures.  From  Equation  (8),  values  of  greater  than  a  decision  threshold  fail  to  reject 
the  hypothesis  that  an  event  is  a  single-point  fully  contained  explosion. 

Statistical  Properties  of  the  Maximum  of  Magnitude  Residuals 

The  energy  radiation  pattern  of  surface  waves  from  a  single-point  fully  contained  explosion  can  be  reasonably 
assumed  to  be  azimuthally  isotropic.  With  f=f  ms  and  cf=cf  /&,  Equation  (6)  gives  the  Ms  covariance  matrix  for  an 
event  observed  by  ny  stations.  Make  the  assumption  that  corrections  are  of  sufficient  high  quality  to  give  f  \ls  £  0. 
With  this  assumption,  the  observed  station  Ms  are  independent  (uncorrelated)  for  an  event,  and  the  Ms  population 
mean  is  zero  for  explosions  and  greater  than  zero  for  earthquakes.  For  source  /,  event  j,  Max{Msip}  <  m  <=>  { Msyk < 
m;  k=  1,  n,j}.  Linder  the  null  (HO)  hypothesis  that  a  seismic  event  is  a  single-point  fully  contained  explosion,  surface 
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wave  energy  is  azimuthally  isotropic  and  observed  station  Ms  have  the  same  probability  distribution  and  so  the  110 
likelihood  is 

P(Msl  '  m . MsUij  <  m  |  //0) 

l\Mst  <m  |  //())•  ••/>(  .)/.<„, .,  <in  |  II 0)  P(Msk  <  in  |  II 0)n  \  (9) 
Order  Statistics  Discriminant  Formulation 

The  energy  radiation  pattern  of  surface  waves  from  an  earthquake  is  not  azimuthally  isotropic.  Some  stations  could 
be  on  an  energy  lobe  and  some  could  lie  on  an  energy  node.  1  n  an  event  identification  setting  the  energy  radiation 
pattern  of  an  event  will  not  be  known.  This  means  that  a  network  average  of  station  Ms  values  could  be  strongly 
influenced  by  nodal  stations.  For  example,  an  earthquake  with  four  nodal  stations  and  one  station  on  an  energy  lobe 
could  look  like  an  explosion  with  Ms  versus  mb  network  average  discrimination. 

The  right  hand  side  of  Equation  (9)  provides  the  likelihood  of  the  random  variable  Max { Ms  ljk f  under  the  null 
hypothesis  of  single-point  fully  contained  explosion.  A  conservative  alternative  hypothesis  (Ha)  likelihood  is 
provided  by 

P(Msk  <  in  I  II a)  x  H(Msk  <  in  |  //0)n,J  1 

that  is,  an  Ms  earthquake  model  at  one  station  and  an  Ms  explosion  model  for  the  other  stations.  The 
Neyman-Pearson  likelihood  ratio  test  reduces  to;  Reject  110  if 

rrvra  I  y'OnXn.j-  —  I)  >  (11) 

where  c,  i| i  and  cp(i/i)  are  a  function  of  the  null  and  alternative  hypothesis  parameter  values.  The  function  tp(ni)  e 
(0,1)  and  monotonic  increasing,  and  so 

«Vm  f  vMK-  I) 


From  Equation  (12),  for  a  given  value  of/:  there  are  values  of  m=Max{Msuk}  that  will  fail  to  reject  110  with  a 
conservative  hypothesis  test  based  on  right  side  of  Equation  (12),  yet  reject  if  the  test  is  based  on  the 
Neyman-Pearson  test  (left  side  of  Equation  (12).  This  is  graphically  shown  in  Figure  2.  The  conservative  hypothesis 
test  is  of  the  form  reject  110  if  Max{Ms-,jk }  >  k,  and  is  conservative  in  the  sense  of  “miss  no  explosions”.  Max{  Msi,n} 
as  a  discriminant  offers  a  cost  efficient  approach  to  discrimination  in  settings  where  path  and  source  corrections 
are  mature. 
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Figure  2.  Graphic  of  the  Ncyman-Pcarson  (black)  and  conservative  (red)  likelihood  ratio  X.  The  blue 
horizontal  decision  line  projects  a  decision  region  for  the  discrimination  Max(Ms).  Values  of 
Max(Ms)  interior  to  this  decision  region  reject  HO  with  Ncyman-Pcarson  test  and  fail  to  reject  HO 
with  the  conservative  test. 


CONCLUSIONS  AND  RECOMMENDATIONS 

The  mathematical  statistics  of  two  new  approaches  to  Ms  versus  mb  discrimination  are  developed  in  this  paper.  Both 
methods  are  designed  to  utilize  advances  in  source  physics  theory  by  using  new  source  models  to  predict  (correct) 
observed  surface  wave  magnitudes.  The  network  average  discriminant  formulation  additionally  includes  source 
model  error  in  the  standard  error  (SE)  of  a  discriminant  effectively  placing  a  lower  bound  on  the  discriminant  SE. 
This  new  property  of  the  Ms  versus  mb  discriminant  correctly  reduces  the  SE  only  through  network  averaging  of 
magnitudes.  The  order  statistics  discriminant  formulation  is  potentially  robust  to  station  coverage  for  an  event  by 
basing  discrimination  on  the  maximum  observed  station  Ms  value.  Validation  analysis  for  these  new  formulations  is 
in  progress  using  global  events  reported  in  reported  in  Bonner  et  al.  (2003),  Bonner  et  al.  (2006),  Marshall  and 
Bashman  ( 1 972).  The  Bonner  et  al.  (2003)  and  Bonner  et  al.  (2006)  papers  research  and  develop  new  signal 
processing  theories  for  Ms  that  are  design  to  Rayleigh  wave  energy  from  small  magnitude  events  (stations  close  to 
source). 
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